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Abstract 



The non-linear hydro dynamic stability of thin, compressible, Keplerian disks 
is studied on the large two-dimensional compressible scale, using a high- 
order accuracy spectral method. We find that purely hydrodynamical per- 
turbations, can develop initially into either sheared disturbances or coher- 
ent vortices. However, the perturbations decay and do not evolve into a 
self-sustained turbulence. Temporarily, due to an inverse cascade of energy, 
which is characteristic of two-dimensional flows, energy is being transferred 
to the largest scale mode before being dissipated. In the case of an inner 
reflecting boundary condition, it is found that the innermost disk is glob- 
ally unstable to non-axisymmetric modes which can evolve into turbulence. 
However, this turbulence cannot play a significant role in angular momentum 
transport in disks. 

Subject Headings: accretion, accretion disks - hydrodynamics - instabili- 
ties - turbulence 
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1 Introduction 



Recently, it has been demonstrated (Balbus & Hawley 1991; see also their 
review Balbus & Hawley 1998) that the origin of turbulence and angular mo- 
mentum transport in accretion disks is very likely related to a linear, local, 
magnetohydrodynamic (MHD) instability-the magnetorotational instability 
(Velikhov 1959; Chandrasekhar 1960; Balbus & Hawley 1991). This instabil- 
ity has been shown to successfully reproduce the order of magnitude of the 
effective viscosity in ionized disks (with a viscosity parameter a 0.001 — 0.1; 
Hawley, Gammie & Balbus 1995, Brandenburg et al. 1995; Vishniac & Bran- 
denburg 1997). The situation in dense and cool disks, like the ones found in 
Young Stellar Objects, is somewhat more ambiguous (e.g. Blaes & Balbus 
1994; Regos 1997) and it has been suggested that the source of turbulence and 
angular momentum in these systems might be due to a hydro dynamical in- 
stability (e.g. Zahn 1990; Dubrulle 1993). However, no linear hydrodynamic 
instability has been identified (e.g. Stewart 1975), and efforts have focused on 
nonlinear instabilities (e.g. Zahn 1990; Dubrulle & Knobloch 1992; Dubrulle 
1992, 1993). The latter authors have derived some interesting properties of 
the turbulence (and the viscosity in the solar nebula), however, they actually 
assumed the turbulence to exist ab initio, rather than rigorously demonstrat- 
ing that it is the result of some instability (basing their assumption on the 
fact that Couette and Poisseuille flows do exhibit nonlinear instabilities in 
laboratory experiments; e.g. Coles 1965; Daviaud, Hegseth & Berge 1992; 
see also Zahn 1990). In fact, until now, it has not been shown that a Ke- 
plerian disk is unstable to a local nonlinear hydrodynamic instability that 
leads to turbulence. Quite to the contrary, using a relatively low-resolution 
finite-difference code, Balbus, Hawley & Stone (1996) have shown that Kep- 
lerian flows are stable to finite amplitude perturbations on the small three- 
dimensional scale, using the shearing-box approximation. They explain their 
numerical findings by pointing out that the epicyclic term is always a sink 
in the energy equation, and consequently, the Reynolds stress is expected to 
decay, even if turbulence transports angular momentum outwards. However, 
no rigorous analytical proof is provided. 

In order to further test the hydrodynamic stability of Keplerian disks, we 
propose in this work to check the non-linear stability of the flow on the large 
compressible two-dimensional scale, since on the small three-dimensional 
scale the flow is expected to be stable (Balbus et al. 1996). The moti- 
vations behind this particular calculation are multiple: (i) two-dimensional 
compressible rotating flows (e.g. the shallow water approximation) do exhibit 
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a transition to turbulence due to the propagation of waves (e.g. Pedlosky 
1987; Cushman-Roisin 1994; see also Clio & Polvani 1996a,b, and see §2.5). 
(ii) The wavelength and the amplitude of the perturbation needed for a sub- 
critical transition to turbulence are not known a priori, (iii) The present 
method of calculation allows us to study the entire disk, thus including the 
effects of the amplification of propagating waves, rotation, curvature and 
compressibility, which might have crucial effects on the flow (see e.g. Cho & 
Polvani 1996a for such a discussion). 

Since the transition to turbulence in shear flows is usually assumed to be 
a strictly three-dimensional effect (e.g. Balbus & Hawley 1998), we discuss 
in the next section the transition to turbulence in three- as well as two- 
dimensional rotating flows. The equations and the numerical method are 
presented in §3, the results are shown in §4, and a discussion follows. 



2 Transition to turbulence 

Planar shear flows are unstable to finite-amplitude perturbations only in 
three dimensions (3D) (e.g. Bayly, Orszag & Herbert 1988), and numerical 
simulations have also shown that a subcritical transition to turbulence in the 
plane Couette flow occurs only in 3D (e.g. Orszag & Kells 1980; Dubrulle 
1991). These results have led to the general assumption that a hydrody- 
namic (subcritical) transition to turbulence in accretion disks can take place 
only in 3D. However, this assumption is based entirely on experiments and 
simulations with incompressible 3D flows, and therefore great caution should 
be exercised in attempts to draw conclusions for 2D compressible flows in 
accretion disks. In particular, waves play a crucial role in the formation 
of turbulence in the shallow water approximation of planetary atmospheres 
(e.g. Cho & Polvani 1996a, 1996b). Therefore, there are good reasons to 
believe that turbulence can in principle at least develop in two-dimensional 
disks as well. 

In this section we review the subcritical transition to turbulence in three- 
dimensional incompressible and two-dimensional compressible flows. 
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2.1 Three-dimensional incompressible flows 

The transition to turbulence in three-dimensional incompressible flows has 
often been considered to be due to an inflectional instability. However, more 
recently it has been realized (e.g. Yang 1987; Hamilton & Abernathy 1994; 
Dauchot & Daviaud 1995; Hegseth 1996) that the destabilization process 
that leads to turbulence is associated with the presence of streamwise vor- 
tices rather than with inflection points: non-transitional flows were found 
to exhibit inflectional profiles as well. A streamwise vortex always induces 
inflections in the streamwise velocity profiles, but the transition to turbu- 
lence occurs only when the strength (circulation) of the vortex is above some 
critical value (Hamilton & Abernathy 1994). The turbulence is fed by the 
streamwise vortex which transfers energy from the bulk motion of the flow 
to the turbulent region (usually a 'spot'). Streamwise vortices in the flow 
are purely three dimensional, and cannot be reproduced by two dimensional 
simulations. 

2.2 Two-dimensional compressible flows 

Two-dimensional flows with inflectional velocity profiles satisfy the necessary 
(but not sufficient) condition for instability (e.g. Fj0rtoft 1950). Therefore, 
one of the features we might expect to observe if the flow is unstable, is 
the presence of inflection points, since this is a necessary (but not sufficient) 
condition for instability in two-dimensional flows. 

However, very little is known about the dynamics of rotating flows when 
compressibility is important, since only a few studies have been carried out 
for supersonic rotating flows in two dimensions (e.g. Farge & Sadourny 
1989; Tomasini, Dolez & Leorat 1996; Godon 1998). In their work on two- 
dimensional transonic shear flows, Tomasini et al. (1996) did not rule out 
turbulence, although it was not obtained numerically. The turbulence found 
in Godon (1998) was the non-linear consequence of the Papaloizou-Pringle 
(1984) instability, which itself was due to an inner reflecting boundary con- 
dition (Gat & Livio 1992). 

In the shallow water approximation of two-dimensional planetary atmo- 
spheres, the shear flow is unstable to wavelengths that are longer than a 
critical value. These waves have a phase speed that matches the mean cur- 
rent velocity within the flow, which permits a transfer of energy from the 
current to the wave. This process is known in geophysical fluid dynamics as 
the barotropic instability (e.g. Pedlosky 1987; Cushman-Roisin 1994). This 
instability is linear, and it is responsible for the formation of turbulence in 
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the atmosphere on the large two-dimensional scales which transport angular 
momentum. An important point to mention here is that two-dimensional 
compressible polytropic flows (with a polytropic index n=l) can be repre- 
sented by the shallow water equations, when the depth h of the shallow water 
flow is replaced by the surface density £, and the surface gravity waves speed 
\fgh corresponds to the sound speed (see e.g. Tomasini et al. 1996; Ingersoll 
& Cho 1997, private communication; Bracco et al. 1998, 1999). Moreover, 
in the thin disk approximation, like in the atmosphere of a planet, the long 
wavelength modes are two-dimensional because of the finite vertical extent 
of the flow and because of the rotation. On the smaller scales (say A < H in 
disks and A < h in the atmosphere) the subsonic flow is not affected by the 
rotation or the curvature (see also Dubrulle & Valdettaro 1992). However, 
while there are many similarities between the Shallow Water Equations and 
the two-dimensional compressible equations for a disk, there are also some 
important differences (see §5). 

3 Numerical modeling 
3.1 The Equations 

The equations (e.g. Tassoul 1978) are written in a cylindrical coordinates 
system (r,(f),z), and are further integrated in the (z) vertical direction (e.g. 
Pringle 1981). We define a mean density p using the definition of the surface 
density E = f* H pdz = J^^ max pdz, where H max has a constant value H max > 
max(H). This relation always holds since by definition the integral vanishes 
for H > H max . The surface density E is then substituted in the equations 
by 2H max p and the equations are further divided by 2H max . The equations 
are written for the density p (where for simplicity we have dropped the bar), 
the radial momentum U = pv r and the angular momentum W = pv^. 
The equations are: the conservation of mass, 
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where Rij (i = r, 0; j = r, 0) is the Reynolds stress tensor. 
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The pressure is given assuming a polytropic relation 

P = Kp l+1/n . 



(4) 



We chose n = 3 for the polytropic index, while the polytropic constant K 
was fixed by choosing H/r at the outer radial boundary. 

We use an alpha prescription (Shakura & Sunyaev 1973) for the viscosity 
law [y = ac s H = acf/fi^), where c s is the sound speed and H = c s /Qk is 
the vertical height of the disk. 

The Reynolds number in the flow is given by 
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where L is a characteristic dimension of the computational domain, v is the 
velocity of the flow (or more precisely the velocity change in the flow over a 
distance L) and v is the viscosity. Since we are solving for the entire disk, 
L pa r and v ~ vk, and the Reynolds number becomes 
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3.2 Numerical method 



We use a Fourier- Cheby she v spectral method of collocation (described in 
Godon 1997) to study waves in accretion disks. Spectral Methods are global 
and of order N/2 in space (where N is the number of grid points; see also 
Gottlieb & Orszag 1977; Voigt, Gottlieb & Hussaini 1984; Canuto et al. 1988) 
and they make use of fast Fourier transforms. These methods are relatively 
fast and accurate and consequently they are frequently used to solve for tur- 
bulent flows (e.g. She, Jackson & Orszag 1991; Cho & Polvani 1996a). All 
the details on the numerical method can be found in Godon (1997), here we 
only give a brief outline of the numerical modeling. We would like however 
to remark that because spectral methods are high-order-accuracy methods, 
the number of points needed to achieve a given accuracy is much less than in 
usual finite difference methods. For example, a modest number of 17 points in 
a spectral code has the accuracy of a low order finite difference solver with a 
100 points (Gottlieb & Orszag 1977; Voigt, Gottlieb & Hussaini 1984; Canuto 
et al. 1988). It is quite remarkable that Orszag & Kells (1980) were able to 
study the subcritical transition to turbulence in the Plane Couette Flow with 
a spectral method with a resolution as low as 16 3 points. Spectral methods 
reach a high accuracy for smooth solutions, and the series do not converge for 
discontinuous solutions (e.g. very strong shocks in the flow can be respon- 
sible for numerical oscillations, the Gibbs phenomenon). However, in the 
present work the shocks that may form are due to the amplification of waves 
propagating inwards in the disk. These shocks are (relatively) weak and they 
do not affect the accuracy of the results presented here. For example, similar 
shocks were obtained due to the amplification of tidally induced spiral waves 
propagating inwards in the disk (with a Mach number of 1.3, Godon 1997). 
These waves were found to reach the inner boundary, even assuming an alpha 
viscosity parameter a = 0.01 and a resolution of 128x32 (Godon 1997), while 
in inviscid models using finite difference schemes (e.g. Savonije, Papaloizou 
& Lin 1994) the tidal waves were damped before reaching the inner boundary 
even though the resolution was higher (200x64). These results show that in 
the circumstances that apply to the present calculations (perturbations of 
a disk) our spectral code has a higher accuracy than usual finite difference 
codes with the same number of grid points. 

A Spectral filter was used to cut-off high frequencies. This implementa- 
tion was used for numerical convenience. It eliminates the high frequencies 
(e.g. due to shocks and aliasing errors of the non-linear terms) which can 
cause numerical instability, while keeping a high enough number of terms in 
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the spectral expansion to resolve the fine structure of the flow (Godon 1997, 
1998; see also Don & Gottlieb 1990). 

In all the numerical simulations, the alpha viscosity prescription was 
introduced in order to dissipate the energy in the small scales (say, for 
wavenumbers k > k u ), and the spectral filter was used to cut-off the highest 
frequencies (for wavenumbers k > ko > k v ) which can grow due to alias- 
ing errors and the Gibbs phenomenon (Gottlieb & Orszag 1977; Voigt et al. 
1984; Canuto et al. 1988). To make sure that the flow is numerically stable 
we check the energy spectra (E k ) of the flow and make sure that the energy 
decreases (say, like Ek oc k" 13 with 3 < j3 < 4) for wavenumbers k > 3M/4, 
where M is the total number of modes (i.e. k = M corresponds to the high- 
est frequency). For low amplitude perturbations, a spectral filter is sufficient 
to damp the high frequencies (by choosing ko = 3M/4). However, when the 
perturbations are strong (like in the present work), a small viscosity is needed 
in addition to the spectral filters (the exponential cut-off filter used here is a 
weak filter). 

4 Results 

We have run several classes of models. In models of the first class, we simply 
perturbed the disk with a finite-amplitude perturbation, and free boundary 
conditions were applied at the inner and outer edges. In the second class of 
models, in addition to perturbing the disk, we also included (throughout the 
run) the tidal effects from a companion star. In the third class we introduced 
in the equations an artificial change in the epicyclic term (see §4.3), and 
finally, in the fourth class, we applied a reflective boundary condition at 
the inner edge of the disk. In all the models the initial rotation law is 
Keplerian and the radial velocity is set to zero. The initial density profile 
was chosen to fit the standard thin disk model (Shakura & Sunyaev, 1973; 
i.e. p(r) oc r~ 15//s ), but models with a flat density profile were also run with 
similar results. 

4.1 Finite perturbations, free boundaries 

In this class of models, we have run more than 10 models, changing the resolu- 
tion (N x Q, where iV is the number of collocation points in the Chebyshev 
expansion for the radial dimension, and Q is the number of points in the 
Fourier expansion for the azimuthal dimension), the ratio of the outer to the 
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inner radius (R ou t/ Rin), the viscosity parameter (a) and the Mach number in 
the flow (M. = (H/r)^ 1 ). The viscosity parameter was varied from a = 0.1 
down to a = 5 x 10~ 5 (more than three orders of magnitude change in the 
Reynolds number) and we tried mainly two values of the Mach number, cor- 
responding to H/r = 0.05 and H/r = 0.15. Not all the computed models are 
presented here. 

The optimal results (in terms of accuracy and speed) were obtained for 
Rout/ Rin = 2 and N x Q = 128 x 128. Given the fact that we are working 
with a high-order accuracy spectral method, the resolution achieved with 
128 x 128 points is equivalent to that obtained by a usual second order accu- 
racy finite-difference method with more than 500 2 grid points (Canuto et al. 
1988). This improved resolution is achieved because the spatial resolution is 
not only a direct function of the number of grid points, but it is also an inverse 
function of the numerical dissipation, which is responsible for the flattening 
and spreading of fine structures in the flow. High order accuracy spectral 
methods have very little numerical dissipation in comparison to standard low 
order accuracy finite-difference methods with the same number of grid points. 

We chose in the present work to perturb the disk velocities with a high 
order mode. The initial perturbation was randomly chosen to have a maxi- 
mum amplitude of about 0.1 — 0.2 of the sound speed. However, even larger 
amplitude perturbations were obtained naturally, due to the amplification of 
the waves as they propagated inward into the disk, resulting in some cases in 
shocks in the inner region of the computational domain. The perturbations 
were placed on a grid as shown in figure 1. The kinetic energy of the initial 
perturbations ranged from 10~ 5 to 10~ 3 of the background kinetic energy of 
the Keplerian flow. 

Figure 1: A grayscale of the density shows the location of the perturbations 
induced in the density of the disk. 

In the first model (model 1) we chose a = 0.1 with R out /R in = 5, 
H/r = 0.05 and a resolution of 64x64. On a very short time scale the 
initial perturbation was damped and, as expected for such a large viscosity, 
the disk reached a stable configuration even before a significant propagation 
of the waves has been noticed. 

In a second model (model 2) we reduced the viscosity parameter to 
a = 0.01 and changed the numerical resolution to N x Q = 128 x 128 with 
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Rout/Rin = 2, keeping H/r = 0.05. In this case the perturbation had time 
to propagate and form spiral waves sheared by the flow. However, the short 
wavelength waves were rapidly damped (see figure 2). After about 5 Keple- 
rian orbits all the waves were smoothed out. We found that in the models 
the waves dissipated on a time scale r dis (r^P^) 1 / 2 , where r v = r 2 jv is the 
viscous time scale and Pk = 2tt/Qk is the Keplerian Period. 

Figure 2: The initial perturbation of the disk induces waves that form spi- 
rals due to the shear in the flow. In this model H/r = 1/20 and the viscosity 
parameter is a = 0.01. The waves are smoothed out and decay on a time scale 
Tdis ~ (tvPr) 1 ^ 2 , where Pk = 2tt/Qk is the Keplerian period and t v = r 2 jv is the 
viscous time scale. 

In model 3 we have further reduced the viscosity to a = 10~ 3 , while keep- 
ing all the other parameters constant. In this case the waves decayed on a 
much longer time scale, and the shear had sufficient time to act strongly on 
the waves. This resulted in the formation of tightly wound spirals in the disk 
(figure 3, the radial wavelengths are very short and are barely visible in the 
inner disk). Here again the waves dissipated after a characteristic time r^ s . 

Figure 3: The same as Fig. 2, with a = 10~ 3 (see text). 

As we have further reduced the viscosity in the disk, the wavelengths of 
the waves in the radial direction decreased to a point where a higher reso- 
lution was needed. At this stage we decided to increase the ratio H/r while 
keeping the same resolution. In model 4 we therefore took H/r = 0.15 and 
a = 5 x 10 -5 . Here the largest waves dissipated on a time scale of the order 
of dozens of orbits, and the effect of the viscosity was to damp the shortest 
wavelengths first. As the longer wavelength modes propagated inward, they 
amplified and formed shocks similar to tidally induced spiral waves. The spi- 
ral waves were not tightly wound and the wavelength in the radial direction 
was not very short due to the lower Mach number in the disk (Figs. 4 and 5). 

Figure 4- The grayscale of the density is shown (at t ps 2 orbits) for model 4 
with a very low viscosity (a = 5 x lO -5 ^ and a lower Mach number (H/r m 1/7). 

Figure 5: A grayscale of the density for model 4 around t w 10 orbits. The dis- 
sipation time of the perturbations is long, of the order of dozens of orbits, however 
the shortest wavelength disturbances disappear within about 10 orbits (see text). 
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As an additional check of the stability of the disk, we have also looked 
at the absolute vorticity of the flow. It is well known from the theory of 2D 
turbulent flows that an inflection point in the vorticity is a necessary condi- 
tion for the appearance of an instability. The absolute vorticity of the flow is 
shown in figure 6, a few orbits after the initial perturbation of the disk. The 
striking feature in figure 6 is a multitude of Z-like shapes criss-crossing the 
whole disk. The Z-shape in the absolute vorticity is a characteristic signature 
of the instability of sheared disturbances (e.g. Haynes 1987). However, while 
this instability is usually driving turbulence in shear flows, in the present 
Keplerian case it eventually stabilizes, as shown in figure 7. This strongly 
suggests that the instability is deriving its energy from the initial perturba- 
tion rather than from the Keplerian motion of the flow itself. 

Figure 6: A grayscale of the absolute vorticity for model 4 after a few orbits. 
The initial perturbations have evolved into Z-shapes, characteristic of the instabil- 
ity of sheared disturbances. 

Figure 7: A grayscale of the absolute vorticity for model 4 after about 10 orbits. 
The flow has stabilized completely. 

In all the preceding models, the size and amplitude of the perturbation 
induced in the density were significant (Figure 5), especially due to the am- 
plification of waves propagating inward (forming shocks). However, the size 
and amplitude of the perturbation induced in the vorticity profile were in- 
deed very small (Figure 6), and the dissipation of the vorticy perturbation 
might be due to the shear and the viscosity (Figure 11 shows that modes 
higher than m 20 are quickly damped for a = 5 x 1CT 5 ). It is important to 
note that the amplitude of the perturbation needed for the development of 
turbulence depends on the Reynolds number (e.g. Dauchot & Daviaud 1994). 
In fact, the critical amplitude for the non-linear growth of disturbances in 
shear flows follows a scaling law which depends on the shape of the initial 
disturbance, as a result of the competition between the non-linear growth 
of the perturbation and the viscous dissipation mechanism (e.g. Dubrulle 
& Nazarenko 1994). Therefore, for a given resolution, a given configuration 
may need a large amplitude perturbation to give rise to an instability. Con- 
sequently, in order to check the effects of a large perturbation in the vorticity 
field, we decided to run an additional model in which the initial vorticity 
perturbation is large in both size and amplitude. The reasons behind this 
vorticity perturbation are multiple. First, two-dimensional flows are unstable 
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to local extrema of the vorticity (inflexion points) that can lead to turbu- 
lence. Second, two-dimensional planetary atmospheric turbulent flows are 
usually perturbed and studied in the vorticy field rather than in the velocity 
field. Lastly, recent numerical simulations using the shallow water equations 
(Bracco et al. 1998, 1999) have shown that vortices can survive in a Keple- 
rian shear flow for many orbits. 

In this model (model 5) we chose H/r — 0.15, a — 1 x 10~ 4 and a 
polytropic index n = 2.5. The initial vorticity perturbation profile is shown 
in Figure 8. 

Figure 8: A grayscale of the initial vorticity perturbation. For clarity, the 
Keplerian velocity background has been subtracted from the angular velocity. The 
maximum amplitude of the velocity field is 5u ~ 0.2c s , where c s is the sound veloc- 
ity. Dark features represent anticyclonic vorticity perturbations and white features 
represent cyclonic vorticity perturbations. 

During the first few orbits all the cyclonic vorticity perturbations were 
sheared by the flow and dissipated, while the anticyclonic vorticity pertur- 
bations became organized into vortices which later merged together to form 
larger vortices. After 15 orbits, only a few large elongated anticyclonic vor- 
tices were left in the flow (Figure 9). 

Figure 9: A grayscale of the vorticity at t w 15 orbits. Anticyclonic vortices 
have developed and merged together to form large elliptical vortices, while the cy- 
clonic vorticity perturbations have dissipated. 

However, all the vortices were eventually sheared and dissipated due to 
the viscosity and the shear in the flow. The perturbations did not develop 
into self-sustained turbulence. In this model, the vortices had typical ex- 
ponential decay timescales of the order of 50 orbits and they can therefore 
survive in the flow for a few hundreds of orbits. Detailed calculations of the 
stability and lifetime of vortices in Keplerian disks (Godon & Livio, 1999a), 
indicate that the exponential decay time of the amplitude A of the vortices is 
inversely proportional to the alpha viscosity parameter. Namely, A oc e~*/ T , 
where r oc a -1 . Therefore, we find that in the range a = 10~ 4 — 10~ 3 , the 
decay time of the vortices is determined by the viscous dissipation that is 
introduced explicitly into the equations (i.e. the Reynolds stress tensor). 

The decreasing kinetic energy of the fluctuations is another indication for 
the stability of the flow, and it was observed in all the models. In figure 10 
we show the decrease of the kinetic energy as a function of time for one of 



13 



the models (model 4). 



It is also interesting to examine the kinetic energy spectrum of the fluc- 
tuations. The energy spectrum is the same in all the models and it is shown 
in Figure 11 for model 4. The slope of the spectrum in the short wavelengths 
(large k) is very close to -3, due to the viscous dissipation of the modes. The 
rest of the spectrum (for the long wavelengths, i.e. small k) is fairly flat, 
in agreement with simulations of two-dimensional compressible turbulence 
(e.g. Farge & Sadourny 1989) and with the turbulence obtained in the inner 
region of the disk (see §4.4). The fluctuations are characterized by an inverse 
cascade of energy typical of two-dimensional flows, and as a consequence, the 
energy accumulates in the lowest (m=l) mode. The spectrum shows clearly 
the dominance of the m=l (wavenumber k — 1, i.e. Log(k) = 0) mode. 
Similar energy spectra were obtained for all the other models in this class. 

Figure 10: The kinetic energy of the fluctuations is shown as a function of 
time. The energy is given on a logarithmic scale in units of the background Ke- 
plerian energy of the flow. The energy decreases steadily with time, with small 
fluctuations appearing towards the end of the run. The same decrease of energy 
was observed for all the models. The model shown here is model 4- 

Figure 11: The spectrum of the kinetic energy of the fluctuations has been 
integrated in space and averaged over the last orbits. The high frequencies are 
damped by the viscosity with a slope approaching -3. The rest of the spectrum 
is fairly flat in agreement with previous results for compressible 2D turbulent ro- 
tating flows. The low order mode m=l has become dominant due to the inverse 
cascade of energy, characteristic of two-dimensional turbulence. The same spec- 
trum of energy was observed for all the models. The model shown here is model 4- 

4.2 Finite perturbations and tidal effects 

As we explained in the Introduction, it is expected that the epicyclic term 
will always act as a sink in the energy equation (see §4.3), thus causing the 
stresses to decay. It is important however to explore whether other effects 
could serve as sources of energy for the fluctuations. In principle at least, 
such sources could offset the effect of the epicyclic term, with a resultant 
growth instead of decay. One such potential source could be tidally induced 
spiral density waves (e.g. Frank, King & Raine 1999). We have therefore 
run a model in which a finite perturbation was applied to a disk containing 
a tidally induced wave. We found that the shock waves did not act as a 
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source, rather, they removed energy from the fluctuations, causing a yet 
faster stabilization. 



4.3 Artificial change of the epicyclic term 

In order to further test the stabilizing effect of the epicyclic term (Balbus 
et al. 1996, eq.2.6b), we introduced the following artificial change of the 
epicyclic term in the angular momentum equations. In the angular momen- 
tum equation (eq.3), the second term on the left hand side can be written 
as: 

l*-r (™) = --rWA (5) 
r dr \ p J r dr 

We introduce a departure from the circular flow using the relation = 
— rQ , and obtain 

1 d (UW\ Id, . ^ d , d , o^pvr 

-— r = -—r{pv rU(t> ) + rQ—{pv r ) + — r 2 fi 6 

r Or \ p J r or Or Or r 

The last term on the right hand side is the epicyclic 'sink' term described in 
Balbus et al. (1996). Since in our simulations the disk is very nearly Kep- 
lerian, we can substitute Q = in the epicyclic term and obtain for the 
epicyclic term |/w r f2. In order to modify this term in the angular momentum 
equation (eq.3) we added on the right hand side a term of the form apv r r, 
where a is a positive constant. The value a = 0.5 corresponds to a precise 
cancellation of the (sink) epicyclic term. We found that while models with 
a < 0.5 were stable, even a slight increase in a above 0.5 (e.g. a = 0.501) 
resulted in a violent instability already in one-dimensional models (used to 
construct the two-dimensional models). The result was similar to the one 
we obtained for Rayleigh unstable rotation laws (specific angular momentum 
not increasing outwards). 

The modification of the epicyclic term has no astrophysical meaning and 
it was introduced only to test the hypothesis that this term is in fact the sink 
term that stabilizes the flow. 



4.4 Reflective inner boundary 

We also performed additional tests of a phenomenon already found by Godon 
(1998). In the latter study, a transition to turbulence was obtained in the 
inner region of the disk, due to the development of a high-order mode of the 
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Papaloizou-Pringle instability (Papaloizou & Pringle 1984) in the non-linear 
regime. The Papaloizou-Pringle instability develops in the inner region of 
the disk when the inner boundary is reflective (see also Gat & Livio 1992). 
In our current simulations, a transition to turbulence was observed when 
the viscosity in the disk was taken to be very low (a ~ 0.001), and the 
order of the unstable mode was very high (m m 20). The turbulent inner 
disk is shown in Figure 12. In this particular model the resolution is 256 x 64. 

In realistic astrophysical configurations, a reflective inner boundary could 
be obtained (i) either at the inner edge of the disk where a sharp drop in the 
density in the boundary layer can reflect waves (e.g. Papaloizou & Stanley 
1986), or (ii) at the outer layers of a white dwarf rotating near break up 
in Cataclysmic Variable systems (e.g. Narayan & Popham 1989), where 
evidence of a boundary layer is missing. 

Figure 12: A grayscale of the density shows the development of turbulence in 
the inner region of a thin Keplerian disk. This non-linear instability is global, due 
to the inner reflective boundary. 



5 Discussion 

Using a high-resolution, high-accuracy spectral method we have performed 
two-dimensional simulations of an entire accretion disk. The Keplerian (com- 
pressible) flow was subjected to a variety of finite-amplitude perturbations. 
Our simulations show that such purely hydrodynamic disturbances do not 
evolve into a self-sustained turbulence. We have shown that, for small per- 
turbations, sheared disturbances develop, but they are unable to tap energy 
from the flow and they subsequently decay. Furthermore, we have shown 
that the presence of a tidal field does not enhance the turbulence but rather 
suppresses it. For large size and amplitude vorticity perturbations, coherent 
anticyclonic vortices form and merge together. However, here too the dis- 
turbance does not tap energy from the flow: the vorticity eventually decays 
and stretches due to the viscosity and shear (respectively). Similar results 
were obtained (Bracco et al. 1998, 1999) using the incompressible shallow 
water equations for the disk. In their models, Bracco et al. were limited to 
a polytropic index n = 1, a constant density profile p, and most importantly 
the value of H/r was not specified (Provenzale 1998, private communica- 
tions). In the present work we provide a more realistic approach and confirm 
the results of Bracco et al. (1998, 1999). The relatively long-lived vortices 
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may have some important consequences for the formation of giant planets in 
protoplanetary disks (Adams & Watkins 1995; Tanga et al. 1996; Godon & 
Livio 1999a) by allowing dust particles to rapidly concentrate in their core. 
However, they do not lead to turbulence and angular momentum transport. 
More details on the stability and lifetime of vortices can be found in Godon 
k Livio (1999a). 

We have also confirmed the fact that the epicyclic term in the energy 
equation, which acts as a sink, is the main cause for the stability of Keplerian 
disks. Our calculations have shown that the transfer of energy to the largest 
scale results in a (transient) dominance of the m=l mode. This could have 
important consequences for short period oscillations of disks (Godon & Livio 
1999b). Finally, we have shown that when an inner reflective boundary is as- 
sumed, the inner disk is globally non-linearly unstable to non-axisymmetric 
modes which can evolve towards turbulence (see also Godon 1998; and argu- 
ments by Brandenburg et al. 1995). While this is an interesting theoretical 
result, we do not think that this instability plays a role in global angular mo- 
mentum transport in disks, since the turbulence is confined to the very inner 
disk. In addition, the relevance of an inner reflective boundary condition has 
been only partially justified in some systems. 

To conclude, while a definitive answer will have to await for full-scale three- 
dimensional calculations, our results strengthen the impression that purely 
hydrodynamical instabilities are probably not the source of angular momen- 
tum transport in accretion disks. MHD turbulence on the other hand not 
only transports angular momentum adequately, but also readily lends itself 
to the a-disk formalism (Balbus & Papaloizou 1999). 
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